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Abstract 

We consider the development of Monte Carlo schemes for molecules with Coulomb in- 
teractions. We generalize the classic algorithms of Bird and Nanbu-Babovsky for rarefied gas 
dynamics to the Coulomb case thanks to the approximation introduced by Bobylev and Nanbu 
[4]. Thus, instead of considering the original Boltzmann collision operator, the schemes are 
constructed through the use of an approximated Boltzmann operator. With the above choice 
larger time steps are possible in simulations; moreover the expensive acceptance-rejection pro- 
cedure for collisions is avoided and every particle collides. Error analysis and comparisons 
with the original Bobylev-Nanbu (BN) scheme are performed. The numerical results show 
agreement with the theoretical convergence rate of the approximated Boltzmann operator and 
the better performance of Bird-type schemes with respect to the original scheme. 

Keywords: Coulomb interactions, plasma physics, Boltzmann equation, Landau equation, 
Monte Carlo methods. 

1 Introduction 

When a gas is far from the thermodynamical equilibrium, the description of the system through 
the fluid equation is not satisfactory and its fundamentals properties depend upon the interactions 
of the particles. Collisional phenomena can be distinguished for long-range interactions and short- 
range interactions. Short-range interactions are typically of rarefied gases and they are described 
through the Boltzmann equation [S], while long-range interactions are normally encountered in 
plasmas and modeled through the Landau- Fokker-Planck equation [3]. Nowadays, numerical sim- 
ulations of plasmas are receiving a great deal of attention both in research and in industry thanks 
to the numerous applications directly connected to these phenomena. In addition, there exist 
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many practical situations in which the so-called Coulomb collisions are fundamental for correctly 
describing the plasma dynamics as for instance in magnetic fusion devices. 

However, while the literature on Monte Carlo schemes for short-range forces is wide (pQ, [2J, [5], 
[7], [5], [T3]) and many efficient methods have been developed for treating these problems, on the 
other hand, how to construct efficient Monte Carlo numerical methods for long-range interactions 
like the Coulomb potential field it is still not clear. The present work is a contribution in this 
direction. We observe that recently, some important results have been achieved by Nanbu and 
Bobylev in (15) and [4] on this subject. They proposed a new Monte Carlo numerical method 
which permits to efficiently simulate Coulomb collisions. The performance of this scheme has been 
studied in details by Caflisch et al in []I5] and [13], in comparisons with the classical scheme for 
simulating Coulomb collisions by Takizuka and Abe [T7]. Monte Carlo methods for plasmas have 
also been developed by Wang et al [20] while an interesting hybrid method has been realized and 
numerically tested by Caflisch et al [8] for accelerating the simulation of Coulomb collisions. 

In the case of Coulomb interactions each particle interacts simultaneously with a large number of 
other particles, thus multiple collisional phenomena are involved adding difficulties to the numerical 
description. However, these multiple interactions can be seen as a number of simultaneous binary 
collisions, each of which gives a small contribution to the relaxation process through a small angle 
scattering between particles. [3j. 

The Landau-Fokkcr-Planck equation is a valid substitute to the Boltzmann operator when 
describing this type of systems. More precisely, the Landau-Fokker-Planck equation can be derived 
as an asymptotic form of the Boltzmann equation in the case of a Coulomb potential field, in which 
the large angle deflection of a charged particle in a multiple Coulomb interaction is considered as a 
series of consecutive weak binary collisions [4j[T0]. However, the classic Boltzmann collision integral 
is still able to describe the interactions, but the typical time between two consecutive collisions 
prohibits construction of efficient explicit schemes, since the resulting time step will be in these 
cases too small and an excessive use of the computational resources is needed. Moreover due to 
the small effect of a single interaction such detailed modeling is unnecessary. 

In the present work, starting from an approximation derived by Bobylev and Nanbu [3] for 
the Boltzmann collision operator, we derive a general class of direct Monte Carlo methods in the 
same spirit of Monte Carlo schemes for rarefied gas dynamic. Moreover, keeping separated the 
discretization of the time derivative and the approximation of the collision operator, we perform a 
series of numerical convergence tests for the approximated collision operator in order to show that 
the effective convergence rate coincides with the one hypothesized in [3] . 

The rest of the paper is organized as follow. In Section 2, we introduce the Boltzmann and 
Fokkcr-Planck equations and their properties. In Section 3, we derive the approximated Boltzmann 
operator and the limiting case in which, for small angle scattering, it converges to an approximated 
operator for the Landau-Fokker-Planck equation. Section 4 concerns the construction of direct 
simulation Monte Carlo methods. Several test problems which show the capabilities of the methods, 
the differences and analogies with the original Bobylev-Nanbu (BN) scheme, and the convergence 
rates are presented in Section 5. Some final considerations are discussed in Section 6. 

2 The Boltzmann and the Landau-Fokker-Planck equations 

Consider the Boltzmann equation 



(1) 



df{x,v,t) 
dt 



+ v ■ V x f(x, v, t) + -V v (a(v)f(x, v, t j) = Q(f, /) 



with the initial condition 



(2) 



f(x,v,t = 0) = fo(x,v), 
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where / = f(x,v,t) is a non negative function describing the time evolution of the distribution of 
particles which move with velocity v £ R 3 at the position x £ il C R 3 at time t > 0. The vector 
a(v) represents the acceleration due to the force acting on particles such as gravity, electric field 
or magnetic field. The bilinear operator Q(f, /) describes the binary collisions between particles 
and is given by 



(3) 



Q(f,f) = [ ( B (\q\, *ff) W)f{<) - f{v)f{v*)]dndv* 
Jm Js 2 V M / 



where S 2 is the unit sphere in R 3 , q — v — v*, n £ S 2 the unit normal. The post collisional velocity 
are computed by 

(4) v' = ^(v + v* + \q\n), < = ^(v + v* - \q\n) 

The collision kernel B(\q\,q ■ n/\q\), which characterizes the detail of the interaction, is defined as 

(5) B(\q\, cos 6) = \q\a(\q\,6), (0 < 6 < tt) 

Here cos8 — q ■ n/\q\ and a(q,9) is the collision cross section at the scattering angle 9, that 
correspond to the number of particles scattered per unit time, per unit of incident flux and per unit 
of solid angle. We introduce also the total scattering cross section and the momentum scattering 
cross section that will be used in the remainder of the paper 

/■7T 

(6) atot(\q\) = 2tt / a(\q\,6)sm0de 

Jo 

(7) v m (\q\) = 2tt f a(\q\,0)sm0(l-cos6)dO 

Jo 

In the case of hard sphere molecules the cross section and the collision kernel takes the form 

(8) <r{q,0)=^, B{\q\,0)=*\v-v.\ 
while in the variable hard sphere case we have 

(9) a{q,0)=C a \v-v*\ a -\ B{\q\,6) = C a \v - v*\ a 

with C a and a positive constants. The case a — is referred as Maxwellian gas while for a = 1 we 
recover the hard sphere model. In the case of Coulomb interactions the Rutherford formula holds 

(10) cj{\q\,6)- " n 



4sin 4 (<2/2) 

where bo = e 2 /(4:neom r \v — v*| 2 ), with e the charge of the particle, eo the vacuum permittivity and 
m r the reduced mass, which corresponds to m/2, if the particles are of the same species, with m 
equal to the mass. Observe that the above formula implies that the scattering cross section tends 
to infinity as the angle 9 tends to zero. In order to obtain finite and meaningful values for the 
total and the momentum cross section it is necessary to introduce a cut-off value for the impact 
parameter. The cut-off value is justified by the shielding effect phenomena, leading to the following 
values for the total cross section and the momentum cross section 



(11) a tot (\q\) = n\. 
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(12) a m Qq\)=47rbllogA 

with X d = (^|f) 1/2 the Debye length and A = 8in(fl JL /2) ■ 

In the case of grazing collisions it is possible to derive from the Boltzmann operator the Landau- 
Fokkcr-Planck operator (see [10] for details) 

(is) Q L (f, f) = |£ / r3 l«K(kl)((M 2 )<% - **) x (A _ /(,)/(«.)*;. 

In the next section we will see how it is possible to construct numerical schemes starting from the 
Boltzmann equation which approximate the Landau operator (fT5|) , 

3 The approximated Boltzmann equation 

From now on, we will focus on the space homogeneous equation without force fields. Once the 
collision term is solved, the solution of the full Boltzmann equation can be recovered by computing 
the transport and the force term through a time splitting. 

Although the divergence of the collision integral has been solved with the cut-off of the scattering 
cross section, the simulation of the Boltzmann equation for Coulomb interactions still represent a 
significant challenge, due to the too high computational cost which is necessary to directly simulate 
the equations with time explicit schemes. In fact rewriting Eq. ([T]) in the space homogenous case 
pointing out the gain and loss term 

(14) % = Q + (f, f) - f(v)L[f](v), L[f](v) = a tot (\q\) f \q\f(v.)dv, 

(15) Q + (fJ)= f f B (\q\,^f) f(v')f«)dndv, 

JwJs* V \q\ J 

it is easy to observe that the large value of the total collision cross section forces the time step to 
be small, thus too many steps become necessary to compute the final solution, yielding this scheme 
useless. In fact discretizing the time derivative we obtain 

(16) f{v, t + At) = AtQ+(f, /) + f(v, t) (l - Ata tot (\q\) J \q\f(v,)dv^ 

now if we want to preserve a probabilistic interpretation we need the coefficients to be positive, 
thus At has to be extremely small if cr tot (|<7|) is very large. 

Recently an approximated Boltzmann operator has been developed by Bobylcv and Nanbu ([!]), 
which permits use of larger time steps during the simulation even in the case of Coulomb collisions. 
Here we try to generalize this approach in order to construct Direct Monte Carlo schemes for small 
particles interactions. 

Rewrite equation (|T|) in the homogenous case in the following form 



(17) %= [ JF(U,q)dv* 

where U = (v + u*)/2 denotes the center of mass velocity, and 

(18) F(U, q) = f(U + q/2)f(U - q/2) = /(«)/(«,) 
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while the operator J is defined as 

(19) JF(U,\q\u)= [ B(\q\,u-n)[F(U,\q\n)-F(U,\q\uj)]dn 

Js 2 

with u) = q/\q\. If we approximate the operator J in equation (JT7J) by 

(20) J= -(exp(Tj) - I) 

r 

where / is the identity operator and t are assumed to be small, the equation reads 

(21) % = - [ (exp(rJ) - T)F(U, q)dv* = - (P T (/, /) - gf) 

at t j R3 t 

with 

(22) PA.f,f)= f exp(rJ)f(v)f(v*)dv* 
The operator exp(rJ) can be written as 

(23) exp(rJ)^(w) = / B T (u} ■ n, \q\)tj)(n)dn 

Js 2 

where i/}(u>) is an arbitrary function and 

00 91 4- 1 

(24) B T (u ■ n, \q\) = £ exp(-A i (|g|)r)P i (w • n) 

is the Green function, with Pi(uj ■ n) the Legendre polynomial and Aj(|g|) equal to 

(25) A,(|«|) = 2tt £ B{p, \q\){l - Pl(j*))dn 
where /i = uj ■ n, — 1 < /x < 1 Using the above expression we obtain 

(26) PAfJ)= f B T (\q\,^)f(v')f(vi)dndv* 



■■xs 2 \Q\ 



Note that 



(27) / S T (w.n,|g|) = l 

Js 2 

3.1 A first order approximation for the Landau- Fokker-Planck equation 

Assume now that the scattering cross section cr(|<7|,#) is concentrated at small angle near 9 w 0, 
thus B T (\q\,^i) is concentrated near /i = 1. In that situation it is possible to derive the following 
formal approximation 

(28) A/(u)~2tt/ B T (\q\,fx)(l-Pi(l) + (l-fi)Pl{l))dfi = 7rl(l + l) [ B T {\q\,n)(l - fi)dfx 
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where P/(l) = 1(1 + l)/2. The approximate Green function reads 

(29) B r (n, \q\) ~ B T V k|) = £ ^ W cxp (-^y^kk m (k|)r) 

The superscript L in equation (129[) means that equation (|2ip with the above kernel approximates 
the Landau- Fokker-Planck equation. For a formal proof we refer to the paper of Nanbu and Bobylev 

(El)- 

Consider now the case of a Coulomb potential field in a single component gas or plasma. This 
choice, with the cut-off of the scattering angle introduced in the previous section, leads to the 
following approximated equation of order 0(r) 

(30) % = -{! D(^,-^-)f(v',t)f(vi,t)dndv*- Q f(v,t)) 



where 



(31) — =4tt. , 

n \4ire m r J {q^ 



1 / e 2 \ 2 gin A 



and 

(32) D{jn, T ) = f T Lip '(A 1 ) *M-IQ + 1)td). 



" 22 + 1 



47T 



4 DSMC schemes for Coulomb Interactions 

Note that is not necessary to work with the collisional kernel D(/j,,tq) computed above, instead 
a simpler function D*(/i, To) can be used, preserving the same accuracy 0(r), if the following 
condition remain satisfied 

(33) £>*(//, to) > 0, 2ir [ D*(fi, r )dfi = 1 



(34) lim D*( l i,T ) = ±5(l-[i) 



(35) lim — / [D*( f i,r o )-D(^T O )}P l ( f i)d^ = 0. 
One possible substitution is represented by 

(36) D *^^ = 4^mhl eXp ^ 
where A — A(t) satisfy 

1 r_ 

(37) coth A — — = exp ^ 



G 



It is now clear that it is possible to apply slightly modified versions of the standard direct Monte 
Carlo algorithms for Maxwell molecules to the equation 



(38) 



with 



(39) 



V*(/,/) = / D* f(v',t)M,t)dndv* 



The only difference is the way the angle is sampled. In most of the DSMC methods (Hard sphere or 
Variable Hard Sphere scattering models) the angle is sampled uniformly over the sphere, while here 
is sampled accordingly to D*(n, To). Let us discretize the time and denote f n (v) the approximation 
of f(v,nAi), the forward Euler scheme can be used to solve Eq. (|38l) 



This equation has the following probabilistic interpretation: a particle with velocity V{ will not 
collide with probability (1 — gAt/r) and it will collide with probability gAt/r accordingly to the 
collision law described by P*(f, /). Observe that the probabilistic interpretation holds till gAt < r, 
otherwise the coefficient in front of f n becomes negative. Note that taking the limit of the above 
relation, gAt = r, leads to the scheme of Nanbu and Bobylev. The possibility to take different 
values of At < t/q permits reduction of the statistical fluctuations and reduction of the error due 
to the time discretization at no additional cost since, in contrast to the Variable Hard sphere case, 
here no acceptance-rejection procedure is present. 

Hence a Monte Carlo algorithm for the solution of the approximated space homogeneous 
Landau-Fokker-Planck equations reads as follows 

Algorithm 1 (Nanbu-Babovsky (NB) for Coulomb Interactions). 

1. Given N samples v® with k = 1,2, ,.,N computed from the initial distribution function 



2. DO n = 1 to utot with utot = t final/ At 
Given {v£,k = 1,..., N} 

(a) Set N c = round(gNAt/2r), where the round is statistical 

(b) Select N c pairs uniformly among all possible pairs 

(c) Perform the collision between i and j particles according to the following collision law 



(40) 




f(v,t = 0) 



i. Compute the cumulative scattering angle cos 9 as 



(41) 



cos e = — ln(exp A +2U sinh A) 



where U is a random number and A = A(t) is computed through the solution of the 
non linear equation 



(42) 



coth A — - = exp 
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ii. With the above value of cos 9 perform the collision between i and j and compute the 
post collisional velocity according to 

(43) v'i = Vi - ^ (q(l - cos 9) + h sin 9) 

(44) v\ = Vj + i(g(l - cos 9) + h sin 9) 

J 2 

where q = Vi — vj , while h is defined as 

h x = q± cos e 

h v = -{lylx cos e + qq z sme)/q± 
h z = - (q z q x cos e - qq y sin e)/q± 
where qj_ = (q y + ql) 1 ^ 2 and e = 2ttU\ with U\ a random number 
Hi. set = V, and v" +1 = v 4 

(d) Set v" +1 = Vi for the particles that have not been selected 
END DO 

We noticed the structure of the approximate operator analyzed in the previous section is similar 
to the structure of the classical Boltzmann collision integral. Thus it is possible to construct, in 
the same spirit of Maxwell molecules for rarefied gas dynamic, a Monte Carlo scheme based on the 
classical Bird method. 

From the inspection of the approximated Landau operator, it follows that the average number 
of significant collisions in a time step At is given by 

(45) 

which means that the average time between two collisions is given by 

The Bird method, in the case of Maxwell molecules, can be seen as a NB scheme in which the 
smallest possible time step At\ — At/N c is used, in fact only one pair collide each Ati 



(47) 



r +1 = (i - ^) r + ^p?(f, f) = (i - 1) r + /) 



Hence the Bird algorithm for the approximated Landau-Fokkcr-Planck equation reads 
Algorithm 2 (Bird for Coulomb Interactions). 

1. Given N samples v® with k = 1,2, ..,7V computed from the initial distribution function 
f(v,t = 0) 

2. set time counter t c = 

3. set At c = 2t/qN 

4. DO n = 1 to n TO r with n TO T = tf ina i/At 
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(a) repeat 

i. Select a random pair uniformly within all possible pairs 

ii. perform the collision accordingly to the collision law defined in the first algorithm 
and produce v i and u ■ 

Hi. Set Vi — v i and Vj = Vj 

iv. update the time counter tc — tc— At c until t c > (n + l)Ai 

(b) Setvl l+1 =Vi,i = l,...,N 

END DO 

The main difference with respect to the previous algorithm is that multiple collisions between 
particles are allowed in At. Moreover in this case the stability condition can be violated and r can 
be greater than gAt. Thus, as the number of samples increase to infinity, the Nanbu-Babovsky 
scheme converge in probability to the discretized approximated Boltzmann equation. On the 
other hand, the Bird scheme converges to the solution of the approximated Boltzmann equation 
increasing the number of samples, in fact At\ approaches zero. Thus for r — > the first converges 
to the solution of the discretized Landau-Fokker-Planck equation while the second converges to 
the exact solution of the Landau-Fokker-Planck equation. 

Remark 1. Observe that t, although in the schemes has a role similar to the Knudsen number in 
rarefied gas dynamics, has not a clear physical meaning. Mathematically the parameter in front of 
the collision operator is a measure of the goodness of the approximation. Once it is fixed it gives a 
model for the interactions between particles which approximates the Landau operator as t goes to 
zero. 



5 Numerical Tests 
5.1 Test case 

The behavior of the DSMC schemes is illustrated through a series of tests in which the relaxation 
of the velocity distribution function with anisotropic temperature T is considered. Thus the initial 
distribution is taken to be ellipsoidal with T x ^ T y — T z . The initial values for temperature and 
density are set to 

(48) T = 4 x Ry, g = 0.5 

where Ry is the Rydberg constant. The initial difference in the temperature is fixed to ATo = 0.8. 
The approximate analytic solution of the Fokker-Planck equation, in the case of small temperature 
difference, for AT = T X — T y is given by [18] 

(49) AT = AT exp" ^ ^ 
the relaxation time tt corresponds to 

(50) — = g e4 lQ g A 

1 ' t t 7 rV^egm 1 /2(fcT)3/2 

where the Coulomb logarithm value is fixed to log A = 0.5. The simulations are run for most of the 
relaxation process; because all the schemes reach the same final equilibrium state and start from the 
same initial data, our interest is to analyze the different behaviors of the methods when particles 
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are both far from these two situations. Thus, fixing tf — 40 with the values above reported we 
obtain ATy/AXb — 0.2 for the analytic solution, in the rest of the equilibrium process the schemes 
and the analytic solution become closer till they coincide. We remark that although analytic, the 
solution is still obtained through approximations and valid for ranges in which ATq is small. This 
is made more clear by the Figures at the end of the section, even for very small r the schemes do 
not relax at the same rate of the approximate formula (|49[) . 



5.2 Simulations 

Our aim is to perform comparisons between the Bird scheme and the Bobylev-Nanbu (BN) scheme. 
The curves which describe the behavior of the Nanbu-Babovsky (NB) scheme with different choices 
of the time step lie between this two extreme cases. In the sequel we will analyze 

• the deterministic error; 

• the statistical fluctuations of the two schemes. 
The sources of errors for the methods are due to 

approximation of the Boltzmann operator; 
finite number of particles; 

• discretization of the time derivative if present; 

• conservative algorithm used for collisions. 

In order to make a fair comparison of the methods and to stress the capacity to describe the 
relaxation phenomena we try to eliminate the common sources of errors. First we compute the 
deterministic error due to the substitution of the original collision operator with its approximation 
and to the discretization of the time derivative. To that aim we increase the number of samples 
and average the solution of M independent realizations removing statistical fluctuations 



(51) u(t) = ^X>(*) 

where u(t) indicates the average solution at time t and 



(52) Ui (t) = 



T(0) x -T(0) y 



with i the realization number. The number of samples used in the convergence analysis test for 
each realization is N=2 x 10 6 while the number of realizations is M — 5. Observe however that, 
using the Bird method together with a large number of particles for each realization, leads to a 
very accurate discretization of the time derivative (the effective time step is a function of 1/N), 
while with the Babovsky-Nanbu scheme the increase of the samples number does not affect the 
treatment of the time derivative. Note that since no acceptance-rejection procedure is necessary 
the two methods have approximately the same computational cost. Summarizing the first test 
tries to measure the deterministic error computing the numerical order of accuracy with respect to 
r of the two methods. From the theoretical analysis we expect both methods to be first order in 
r. Note however that for BN method the error is due to both the approximation of the operator 
and the to discretization of the time derivative. The order of accuracy r in r is computed as 

/ , , x \uUt) -u(2t)\ 
(53) r(r) = log a fl(r), R(r) = ^ - ±- ]' 
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with R the error ratio. Our second purpose is to measure the stochastic fluctuations of the two 
methods. To this aim we compare the two variances defined as 



Fixing the parameter r, i.e. the approximation of the Boltzmann operator, the variances of the 
two methods are compared for increasing number of samples starting from N — 100 to TV = 3200. 
In this test the number of realizations is chosen equal to M — 1500. Thanks to multiple collisions 
we expect the Bird scheme to have slightly less fluctuations with respect to the BN scheme. 

5.3 Results 

Here we report the solution of the tests described in the previous Section. In Figure [T] the solution 
for the relaxation of temperature in the different directions is showed for the Bird method, while 
in Figure [3] for the BN method. The behavior of the schemes has been analyzed using six different 
values for the parameter r 

(55) t = 2t = 1t = 0.5t = 0.25 r = 0.125 r = 0.0625 

Moreover the solution with t = 0.03125 with N — 2 x 10 6 and M — 10 realizations has been 
computed as a reference solution. In both Figures the analytic and approximated solution is 
reported (blue line) showing a discrepancy with the computed solution even with the more accurate 
one. The time step used in the BN scheme is chosen in order to satisfy the relation gAt/r — 1. In 
Figures [2] and 0] the convergence rate r is plotted for respectively Bird and BN. Both the schemes 
approach the value 1 when r — > as expected from the theory. Anyway it is possible to observe 
from Figure [5l in which the solution of the two methods for the same values of r (respectively 
r = 2, r = 1, r = 0.5 and r = 0.25) has been compared, that the two algorithm furnish a different 
relaxation rate for large values of the approximation parameter of the collision operator, while 
for small values the two methods in practice coincide. This behavior can be explained observing 
that while r — > also the time step At — > thus the error introduced by the BN scheme in the 
discretization of the time derivative disappears. 

In order to show the different performance in terms of statistical fluctuations we fixed the 
parameter r = 1 and perform several simulations increasing the number of samples N. In Figure 
[5] a comparison between the two variances, obtained with Bird and BN, has been plotted. The 
statistical fluctuations of the two method are approximately the same for all the initial choices of 
TV; nevertheless it is possible to see how the Bird scheme oscillates slightly less then BN scheme 
for all the values. This behavior is mainly due to the presence of the time discretization error in 
the BN scheme. If, instead of choosing a large value for r we keep it small, the variance of the 
two methods become practically the same. In fact, if we want a fine approximation of the collision 
operator with the Bobylev-Nanbu method the time step has to be very small, which means we are 
neglecting the time discretization error. 

6 Conclusion 

In this work we have proposed a generalization of the Monte Carlo scheme proposed by Bobylev 
and Nanbu for the solution of plasma physics problems in which the predominant collisions are 
of Coulomb type. This result is achieved by extending the classic Nanbu and Bird algorithms for 
rarefied gas dynamic to the case of plasma physics. The new methods provide more accurate results 
for a fixed r with respect to the original BN algorithm without increasing the computational cost. 



(54) 




i=l 
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Relaxation of anisotropic data: (n ,T ,dT /T )-(0.5,4Ry,0.8) Bird Scheme 




t (N,realization)=(2000000,5) 

Figure 1: Relaxation of the velocity distribution function with anisotropic initial data for different 
values of r. Bird scheme. 



In the resulting algorithms not all particles collide at each time step and some particles collide 
more than once. From the physical point of view this result is counterintuitive, in fact, all other 
computational models in literature about Coulomb interactions are based on all particles colliding 
simultaneously. In the limit of small values of t all methods become essentially equivalent. 

In future we hope to extend the methods described in |16) for rarefied gas dynamics to Coulomb 
interactions and to generalize the hybrid techniques developed in [TTJ [T^] to plasma physics prob- 
lems close to thermodynamic equilibrium. Another interesting research direction consists in de- 
veloping a more accurate approximation in r of the Landau operator. This would allow the use of 
larger time steps in the simulations. 
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pointwise convergence rate Bird Scheme 
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Figure 2: Pointwise order of accuracy r(r) = log2(i?(r)) for the Bird scheme for decreasing values 
of the parameter r. 



Relaxation of anisotropic data: (n ,T ,dT /T )-(0.5,4Ry,0.8) Nanbu Scheme 
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Figure 3: Relaxation of the velocity distribution function with anisotropic initial data for different 
values of r. Bobylev-Nanbu scheme. 
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Figure 4: Pointwise order of accuracy r(r) = log2(i?(r)) for the Bobylev-Nanbu scheme for de- 
creasing values of the parameter r. 
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Figure 6: Comparison of the Bird and Bobylev-Nanbu variance for r = 1 and N = 100 top left 
N = 200 top right N = 400 middle left N = 800 middle right N = 1600 bottom left N = 3200 
bottom right with M = 1500 realizations. 



16 



[15] K. Nanbu. Theory of cumulative small-angle collision in plasmas. Physical Review E, 
55(4):4642-4652, 1997. 

[16] L. Pareschi and G. Russo. Time relaxed Monte Carlo methods for the Boltzmann equation. 
SIAM J. Set. Corn-put., 23:1253-1273, 2001. 

[17] T. Takizuka and H. Abe. A binary collision model for plasma simulation with a particle code. 
J. Comp. Phys., 25:205-219, 1977. 

[18] B. A. Trubnikov. Review of Plasma Physics. Consultant Bureau, 1965. 

[19] C. Wang, T. Lin, R. E. Caflisch, B. Cohen, and A. Dimits. Particle simulation of Coulomb 
collisions: Comparing the methods of Takizuka & Abe and Nanbu. J. Comp. Phys., 227:4308- 
4329, 2008. 

[20] W. X. Wang, M. Okamoto, N. Nakajima, and S. Murakami. Vector implementation of non- 
linear Monte Carlo Coulomb collisions. J. Comp. Phys., 128:209-222, 1996. 



17 



